NASA / CR-1999-209 1 20 / VOL4 



Acoustic Treatment Design Scaling Methods 

Volume 4: Numerical Simulation of the Nonlinear 
Acoustic Impedance of a Perforated Plate Single-Degree- 
of-Freedom Resonator Using a Time-Domain Finite 
Difference Method 

R. E. Kraft 

General Electric Aircraft Engines , Cincinnati , Ohio 


National Aeronautics and 
Space Administration 

Prepared for Langley Research Center 
under Contract NAS3-26617, Task 25 


Langley Research Center 
Hampton, Virginia 23681-2199 


April 1999 




Available from: 


NASA Center for AeroSpace Information (CASI) 
7121 Standard Drive 
Hanover, MD 21076-1320 
(301) 621-0390 


National Technical Information Service (NTIS) 
5285 Port Royal Road 
Springfield, VA 22161-2171 
(703) 605-6000 


Abstract 


Single-degree-of-ffeedom resonators consisting of honeycomb cells covered by 
perforated facesheets are widely used as acoustic noise suppression liners in aircraft engine 
ducts. The acoustic resistance and mass reactance of such liners is known to vary with the 
intensity of the sound incident upon the panel. Since the pressure drop across a perforated 
liner facesheet increases quadratically with the flow velocity through the facesheet, this is 
known as the nonlinear resistance effect. In the past, two different empirical frequency 
domain models have been used to predict the Sound Pressure Level effect of the incident 
wave on the perforated liner impedance, one that uses the incident particle velocity in 
isolated narrowbands, and one that models the particle velocity as the overall velocity. In 
the absence of grazing flow, neither frequency domain model is entirely accurate in 
predicting the nonlinear effect that is measured for typical perforated sheets. The time 
domain model is developed in an attempt to understand and improve the model for the 
effect of spectral shape and amplitude of multi-frequency incident sound pressure on the 
liner impedance. This is of particular concern with regard to frequency-scaling of acoustic 
treatment, since the possibility exists that the nonlinear resistance and mass reactance 
effects may be significant compared to the grazing flow effects at high frequencies. A 
computer code for the time-domain finite difference model is developed and predictions 
using the models are compared to current frequency-domain models. 
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Numerical Simulation of the Nonlinear Acoustic Impedance of 
a Perforated Plate Single-Degree-of-Freedom Resonator 
Using a Time-Domain Finite Difference Method 


1. General Theory 

1.1. Introduction 

1.1.1. Background 

Single-degree-of-freedom (SDOF) resonators consisting of thin-walled honeycomb tubular 
cells covered by a porous faceplate are the simplest form of acoustic lining in use as noise 
suppression treatment in aircraft engine ducts. Figure 1 shows the typical construction of an 
SDOF treatment panel, for which the significant parameters are the cavity depth, L, the faceplate 
thickness, t, the orifice hole diameter, d, and the faceplate open area ratio, or porosity, a. The 
design of SDOF liners consists of determining the particular combination of parameters that 
provides the best suppression performance under the design constraints. 



The design is complicated by the fact that SDOF liners covered by perforated orifice 
sheets behave nonlinearly, that is, their acoustic impedance is a function of the Sound Pressure 
Level (SPL) of the sound wave incident on the surface. This nonlinear behavior arises due to the 
resistance to flow through an orifice, which is such that the pressure drop across the orifice is a 
quadratic function of the DC flow velocity through the orifice: 
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( 1 ) 


Ap = AU + BU 2 

where A and B are constants that will be specified later and U is the through-flow velocity 
incident on the facesheet. The acoustic resistance of the resonator, which is the complex ratio of 
incident acoustic pressure to acoustic velocity at the surface of the liner, thus depends on the 
acoustic velocity, which in turn depends on the response of the liner, which is initially unknown. 
An iteration must be performed on the acoustic velocity to determine the acoustic impedance. 
Fortunately, this iteration generally converges quite rapidly. 

The acoustic impedance of SDOF perforated plate resonators is usually predicted using 
frequency-domain models in closed form, with the exception of the velocity iteration 1 . It has been 
found that the impedance at a given frequency depends not only on the acoustic velocity in a 
narrowband frequency range surrounding the frequency, but also on contributions from acoustic 
velocities at other frequencies. Two empirical models currently exist to determine the acoustic 
velocity to be used at each frequency, and the two models are at opposite extremes. The first 
model uses the acoustic velocity in a narrowband around the frequency of calculation, so that the 
acoustic SPL applied is the pressure attributed to that narrowband range. The other model uses 
the square root of the overall sum of acoustic velocities for the entire pressure spectral range, 



where n is the index for the frequency of calculation and N^, is the total number of frequencies. 
The narrowband velocity assumption tends to underpredict the resistance, but the overall velocity 
assumption, which is the more accurate, still tends generally to obscure the finer variations in 
resistance due to pressure spectrum shape. 

Although the frequency-domain models account for the resistance nonlinearity and are 
efficient to apply, they omit an important part of the physical process. The nonlinearity in 
damping causes the acoustic velocity at the faceplate to become non-sinusoidal, thus scattering 
energy into higher harmonics of frequency. This energy scattered from lower to higher 
frequencies will have an effect on the impedance at higher frequencies. To properly account for 
this energy scatter requires a time-domain model. This may be particularly important at very high 
frequencies for scale model treatment panels. 


1.1.2. Objectives 

The objective of this study is to develop a numerical time-domain model that simulates the 
nonlinear impedance effects of an SDOF resonator consisting of a perforated plate over a finite 
depth cavity. Since acoustic impedance is fundamentally a frequency-do main quantity, the model 
converts the acoustic pressure and velocity solution from the time do main to the frequency 
domain to predict the acoustic impedance over a range of frequencies, given an input SPL 
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spectrum. The purpose of this report is to document the analysis, provide a guide to use of the 
computer program, and make a limited preliminary model assessment by comparing time-domain 
predictions, frequency-domain predictions, and measured data. 

The computer impedance prediction model in its current form must be considered to be a 
preliminary version. For increased utility, it will be possible to include more advanced models of 
the velocity or displacement dependence of the orifice end effects, the effects of grazing flow on 
the nonlinear impedance, and the time-dependent effects on the linear resistance term. The 
program will allow computer experiments with various assumptions regarding the physical 
parameters. In particular, it is hoped that the model will allow investigation of effects of 
frequency scaling, that is, behavior at low Reynolds’s numbers and high frequencies. 


1.1.3. Historical Background 

Several researchers have considered the solution to the acoustic properties of a perforated 
plate from a numerical integration point of view. Rice 2 was one of the first to analyze the 
nonlinear behavior of a perforated sheet liner using a time-domain numerical integration approach. 
In his paper, Rice develops a nonlinear resistance model for the effect of multiple frequency 
excitation. The method also includes the effect of steady grazing flow and a high frequency, small 
displacement effect. Rice remarks on the difficulty of including an accurate model for the linear 
resistance term, which is frequency-dependent, and develops an approximate formulation, but 
notes that generally the contribution of the linear resistance term is small in the presence of 
grazing flow for conventional perforates. Rice models a finite depth cavity using a four-element 
lumped mass reactance and stiffiiess reactance model, which allows the cavity response to be 
approximated using transmission line theory rather than using an actual finite difference solution. 

Cummings 3 developed a one-dimensional numerical time-domain model to simulate the 
nonlinear resistance of a perforated plate. Cummings’ models allow several different termination 
conditions on the right side of the faceplate: an infinite tube (pc termination); an infinite baffle 
termination; or a small volume Helmholtz resonator termination (maximum cavity dimension 
much less than a wavelength). These choices of geometry reduce the (nonlinear) integration to an 
integration in time only, eliminating the need for a spatial numerical integration in the spatial 
region to the right of the faceplate. 

Cummings was interested primarily in determining the acoustic transmission loss through 
the perforated plate, as opposed to the impedance at the facesheet. He considers primarily 
transient pressure excitation (tone bursts and pulses). His numerical results show excellent 
agreement with experimental measurements made as part of the study. Cummings also includes 
an empirical model for time variation of the end effect. 

Zorumski and Parrott 4 , using a Fourier decomposition approach, derive a formula for the 
nonlinear effect of high amplitude, multiple-frequency pressure excitation on acoustic impedance 
of a perforate. The method is based on the DC flow resistance measured for the perforate and 
indicates that the impedance at any frequency is a function of the overall acoustic velocity 
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summed over all frequencies. Zorumski and Tester 5 and Rice 6 argue that the particle velocity at 
each frequency should be replaced by the total effective particle velocity summed over all 
frequencies. 

Hersh and Rogers 7 , Hersh and Walker 8 , and Hersh 9 have developed a fluid mechanical 
model for nonlinear orifice behavior in which they assume the particle velocity approaches the 
orifice in a spherical manner. The analysis is limited to a calculation at the t uning frequency of a 
Helmholtz resonator, which is assumed to be known. The analysis is semi-empirical in that it 
requires measured data to determine two key parameters, an effective orifice radius and an 
effective orifice discharge coefficient. The solution for the nonlinear resistance is obtained 
numerically by integrating a second-order differential equation in the time domain. Other than the 
differences in the equations due to the initial orifice flow assumptions, the method of Hersh is 
quite similar to that of Cummings. 

The development of the model in this report is based directly on the work of Rice and on 
the work of Cummings (which is, in turn, based on the preceding work of C umming s and 
Eversman 10 ). The unique element of this effort is the incorporation of a central difference 
numerical integration scheme to account accurately for the effects of a resonator cavity of 
arbitrary depth. 


1.1.4. Assumptions 

The model is one-dimensional in space, so that only a plane wave is assumed to be 
incident upon and propagate within the cavity. The input acoustic pressure on the perforated 
plate surface is given. The acoustic resistance and mass reactance of the faceplate are assumed to 
be lumped in a thin sheet, where the sheet thickness and orifice hole diameter are much less than 
the shortest wavelength considered. 

For this preliminary development, the mass reactance of the perforate orifices is assumed 
to be constant (i.e., the end effects are assumed not to be a function of velocity), and the orifice 
discharge coefficient is assumed to be constant in time. The end effects are incorporated through 
the input of an arbitrary equivalent orifice length. It is assumed that wave attenuation due to 
viscous effects for propagation inside the cavity can be neglected. 


1.2. Description of Acoustic Model 
1.2.1. Model Geometry 

The acoustic model of the SDOF resonator, shown in Figure 2 , is that of a resonator cell 
mounted at the end of a long tube. The diameter of the tube is sufficiently small tha t only the 
plane wave is cut-on (propagates unattenuated) at the highest frequency of interest. The tube is 
assumed to be terminated with a piston source some arbitrary distance from the facesheet on the 
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left. A pressure wave traveling to the right with a known frequency and amplitude content 
impinges on the treatment surface. 



The acoustic properties of the perforated plate are lumped at the plane of the faceplate, 
and a hardwall reflecting surface is located a known depth behind the faceplate. The input to the 
model is the acoustic pressure on the left side of the faceplate, which is combined with the lumped 
faceplate resistance and mass reactance to provide the left boundary condition for the cavity. The 
acoustic one-dimensional wave equation determines the propagation inside the cavity, and a 
hardwall, zero particle displacement condition is the right-hand boundary condition. 


1.3. Acoustic Wave Equation for Particle Displacement 

The partial differential equation for the acoustic wave propagation is written in terms of 
acoustic particle displacement. The acoustic particle velocity, £,, can be related to the “strain” of 
the medium as 


Strain = — (3) 

dx 

where x is the spatial variable. The one-dimensional constitutive equation for the medium can be 
written as 


P = 



where p is acoustic pressure and A is the bulk modulus, which for air is 

A = pc 2 


( 4 ) 


( 5 ) 


5 




where p is the density of the medium and c is the speed of sound. The acoustic particle velocity is 
given by 


v = 


% 

at 


( 6 ) 


The wave equation for acoustic particle displacement inside the cavity is 

a^ = jL^i 

dx 2 c 2 at 2 

This will be expressed in central difference form below. 


1.4. Boundary Condition for Perforated Sheet and Cavity 

At x = L the boundary condition for particle displacement is sim ply 


S = 0 


( 8 ) 


which reflects the hardwall boundary requirement of zero particle displacement (and velocity). 

The acoustic boundary condition at x = 0 is derived from the momentum balance equation 
and the fact that the acoustic velocity is continuous across the facesheet (i.e., acoustic velocity 
must be equal by continuity just to the left and to the right of the facesheet). The pressure applied 
to the fluid in the cavity just to the right of the facesheet will be the incident pressure at the left of 
the facesheet minus the pressure drop across the facesheet due to resistive and inertial effects. 
The pressure to the right of the facesheet must be balanced by the restoring force of the elastic 
medium in the cavity. This can be expressed as: 




at" 


at 


dx. 


( 9 ) 


where M is the effective mass per unit area, and D is the effective damping per unit area. The 
particle displacement £ is the value averaged over the area of the facesheet. 

The effective mass per unit area can be written as 


M = 


Pi 

aC d 


( 10 ) 


where £ is the effective orifice length (sheet thickness plus end effects), a is the open area ratio or 
porosity, defined as the area of the holes divided by the corresponding area of the facesheet, and 
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Cd is the orifice discharge coefficient. The orifice discharge coefficient is defined as the area of 
the vena contracta of the flow from the orifice divided by the geometric area of the orifice. To 
begin, we will consider only the linear part of the facesheet resistance Rl, which can be taken as 
the DC flow coefficient A, from Equation (1), so we take 

D = R l (11) 

The damping force term will later be extended to include the effects of nonlinear terms. 


1.5. Time Scale Transformation and Numerical Stability 

We now make a change to the time scale of the differential equation and boundary 
conditions. If we let 


x = ct 


( 12 ) 


then the wave equation from Equation (7) is transformed to 

d 2 g, = 

dx 2 dx 2 


(13) 


and the boundary condition, Equation (9) becomes 


pfc 2 

aC D 


4_ RL 3 +P(t)+A #=o 


dt 


dt 


dx 


(14) 


We substitute for A from Equation (5) and normalize by pc 2 to get the non-dimensionalized form 
of the boundary condition at x = 0: 


t l a 2 ^ 

^Cd _ dt 2 


Rl^ 

pc dt 


pc 


-p (,) +£=° 

OX 


( 15 ) 


This will be the form of the partial differential equation and boundary conditions to which 
the central difference formulation is applied for the numerical integration. We define a numerical 
integration step size in x as Ax and a step size in t as At. 

As described in numerous references on numerical integration of partial differential 
equations 11 ’ 12,13 ’ 14 , the convergence and stability of the numerical integration is obtained when the 
ratio of Ax to Ax is less than or equal to 1.0, or, choosing the equality condition, then 


Ax = Ax = cAt 


(16) 
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This is equivalent to integrating along a characteristic line of the acoustic wave propagation, and 
is the condition used in the central difference scheme for this analysis. 


1.6. Calculation Scheme and Output 

The numerical integration problem for the SDOF resonator is of the initial value type. The 
medium in the cavity is at rest at time zero, at which time a forcing function (the incident plane 
wave pressure) is applied to the outside surface of the perforated faceplate. An initial transient 
will occur in the cavity as the pressure wave propagates along the cavity and is reflected from the 
hard end. Eventually, a steady-state pressure field response will be achieved in the cavity 
(transient responses will be damped by the facesheet). 

The numerical integration is modeled in space-time as shown in Figure 3, which illustrates 
the x-t computation point grid. The space points are denoted by the index i and the time points by 
the index j, such that the i* spatial position is given by 

Xj=iAx i = 0,1,2,. ..,N (17) 


and the j 01 time value is given by 


Tj=jAt j = 0,1,2,... ,oo ( 18 ) 

Note that there are N+l space divisions in the cavity and that the hardwall is located at 
x = L = NAx. 

The general outline of the computational simulation of SDOF resonator impedance is 
shown in Figure 4. The incident pressure is applied at the outer surface of the perforated sheet by 
“ramping up” the pressure to its final amplitude over a given number of periods of the lowest 
frequency in the pressure signal. The computation is then given a number of longest periods plus 
several cavity traverse times as settling time for the transients. 
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Figure 3 Computational grid for perforated plate/resonator cavity impedance time- 
domain numerical integration model. 


Overview of Numerical Integration Computation of 
Perforated Plate/Resonator Cavity Impedance 



Figure 4. Flowchart for numerical integration computation of perforated 
plate/resonator cavity impedance. 
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After the solution has achieved steady state, the first data acquisition sequence is initiated, 
3 ud the time, pressure and acoustic particle displacement are stored in a block as if they were 
being measured as a digitized time series. The acquisition of data blocks continues until the 
number of blocks specified for frequency averaging has been acquired. The numerical integration 
time step may be much finer than the data acquisition time step, so that there may be many more 
points calculated than acquired for the Fourier analysis. 

Fast Fourier Transforms are performed on the time history data blocks, giving frequency- 
averaged acoustic pressure and particle displacement spectra. The acoustic velocity is computed 
as i© times the acoustic particle displacement, where © = 27tf, and the impedance spectrum is 
found by dividing acoustic pressure by acoustic velocity in each frequency band. 

Output files are obtained for the time histories of the input pressure and resulting particle 
displacement and velocity, as well as for the pressure, acoustic velocity, and impedance spectra. 
The output files can be used as input to standard plotting routines (or spreadsheets) to obtain 
graphical output. 
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2 . 


Finite Difference Equation Formulation - Linear Form 


2.1. Central Difference Formulation for the Partial Differential Equation 

Using the index i to denote x grid points and the index j to denote x grid points, the 
acoustic particle displacement at position i and time j can be written as £,;j. The finite difference 
forms for the first and second derivatives are given by 


dt, _ £i,j+i-£i,j-i 

(19) 

dx 2 Ax 

d £ _ Si+lJ “ Si-Ij 

(20) 

dx 2Ax 

d \ _ Si,j+1 ~ 2 £ij +£i,j-i 
3x 2 (Ax) 2 

(21) 

n 2 / a __\2 

(22) 


dx 2 (Ax) 2 


The central-difference formulation for the partial differential equation (13) can be written 

5i.i,j - 2 4i,j + £i-i,j ey+i-gy-Hkj-i 
(Ax) 2 (Ax) 2 

Using the stability relation from Equation (16), this can be simplified to 

£i,j+l = £i-l,j +£i+l,j ~£i,j-l 

Figure 5 shows graphically how the particle displacement at time j+1 is related to the 
displacement at prior times and locations. 


(23) 


(24) 
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Figure 5. Particle displacement time-step prior influence graph. 


2.2. Central Difference Formulation for the Boundary Conditions 

The finite difference form for the boundary condition at x = L, i = N is sim ply 

,j = 0 

for all j-values, since the acoustic particle displacement is zero at the hard wall. 

The finite difference form of the boundary condition at x = 0, Equation (15), is 
f m R\. 2m (m R\ 1 \ , \ 

U 2+ 2hj 5oj+i h 2 ^ + w » r- H + ^7 Pi + i£ ( ^ - 0 

where we define 


and 


h = At 

i 

m = 

aC D 



pc 


(25) 


(26) 


(27) 

(28) 


( 29 ) 
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Equation (26) can be further simplified to the form 


Cl£o,j+l + ^2^0,j-l - C 3 £oj +C 4 (^l,j _ ^-l,j) + C 6 Pj - 0 


(30) 


where 


C, 


m R 
h^ + 2h 


(31) 


and 


C 2 


m___R 
h 2 2h 


C 3 


2m 


C 4 


2h 



(32) 

(33) 

(34) 


(35) 


(The coefficient C5 has been reserved for later use.) Note that the boundary condition equation 
requires the definition of a fictitious particle displacement at the i = -1 grid point in x, one step to 
the left of the boundary condition surface. 


2.3. Initial Conditions 

The initial conditions on particle displacement are defined along the first horizontal row of 
finite difference grid points, defined by the line j = 0 (for x = 0). Since the medium starts at rest, 
we also require that the acoustic velocity is zero along this line. These conditions are given in 
finite difference form as 


and 


Si,o=<> 

5u-U-i = ° 


(36) 


(37) 
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for all i values. Note that the velocity initial condition requires the specification of the particle 
displacement at the j = -1 time step. 


2.4. Step-Up and Step-Across Algorithms 

To visualize how the numerical integration will proceed, we define two integration 
algorithms, the Step-Across algorithm and the Step-Up algorithm, “across” referring to the space 
dimension and “up” referring to time. 

Consider E, for any succeeding time step at x = 0 (i = 0). Equations (24) and (30) can be 
combined to give 


£o,j+i= - Djpj D 2 £o,j-i + I ) 3^o,j (38) 

where 


Di = 


Cj C 4 

q-c, 

c. 


D-j = 


Cl c 4 


and 


D 4 = 


2C 4 
C] -C 4 


(39) 

(40) 

(41) 

(42) 


This is the Step-Up algorithm, which gives the next time increment from prior time values at the 
faceplate. The Step-Up algorithm, since it is defined at x = 0, incorporates the perforated plate 
boundary condition into the solution. 

With the next highest time value (j+1) given at the left edge, we can then step across to fill 
in the row j+1 (i — 1 to N) using Equation (24), which becomes the Step-Across algorithm. Once 
the first two rows are determined using the initial conditions, the rest of the solution can be found 
for all times. The Step-Across algorithm incorporates the resonator cavity response into the 
solution. 


When the Step-Across algorithm is used to fill in all rows for j = 1 and higher, the zero 
acoustic displacement boundary condition at x = L (i = N) must be enforced. 
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2.5. Starting the Numerical Integration 


The values for the first row are already given in Equation (36). To obtain the second 
row, we use Equation (38) at j = 0 and the initial condition on acoustic velocity in Equation (37) 
at i = 0 to obtain the particle displacement at i = 0 and j = 1, 


D, 


^• ,= ud: Po 


(43) 


The Step- Across algorithm can be used to fill in the rest of row two. Making the 
substitutions, it can be shown that £i,i is zero for all i. This is not surprising, for the acoustic 
perturbation which has just initiated due to turning on the incident pressure at the surface has had 
no time to propagate into the cavity. 
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3. Incorporation of the Nonlinearity in the Boundary Condition 
3.1. Boundary Condition Equation 

We now add a nonlinear resistance terra to the perforated plate pressure-balance boundary 
condition, so that it takes the form 


pc 


-P- 


r l <%> 

R NL 

dZ 

1 

pc dx 

p dx 

dx 

ctC d 


— ^+ — = 0 
dx 2 dx 


(44) 


where Rnl is the nonlinear resistance coefficient which can be identified with the coefficient B in 
the DC flow resistance given in Equation (1). The absolute value of the velocity must be taken to 
assure that positive damping occurs on both sides of the acoustic velocity cycle. 

It is useful to review the physical interpretation of each of the terms in the perforated plate 
boundary condition equation. The first term is due to the applied acoustic pressure, p, on the 
perforate outside surface. The second term is the linear contribution to the damping force, 
proportional to the particle velocity. The third term is the nonlinear contribution to the damping 
force, proportional to the square of the acoustic velocity. Together, the two dampi ng terms are 
equivalent to those that would be obtained from a DC flow resistance measurement, makin g the 
identification of the rms value of the acoustic velocity with the magnitude of the DC flow 
velocity. The fourth term is the inertance contribution of the mass reactance, dependent upon an 
equivalent orifice length, t. The last term is the restoring force due to the compression of the 
elastic medium in the cavity at x = 0. 


3.2. Finite Difference Formulation 


Equation (44) can be expressed in finite difference form as 

£o, J+ l =-D iPj -D 2 ^o,j-l +D 3^0,j +D 5(£o,j+l ~£o,j-l)|£oj+l - £o,j-l 

where Di through D 4 have the same definitions as given previously, 

C5 


(45) 


D 5 = 


c,-c. 


(46) 


and where 


n _ r nl 


(47) 


16 



Ci through C 4 have the same definitions as before. Equation (45) cannot be solved explicitly for 
4 oj+i, and must be solved by an iterative technique. 


3.3. Nonlinear Velocity Iteration Algorithm 

The particle displacement at time j+1 is extracted from the nonlinear equation (45) using a 
Newton-Raphson iteration technique. We write the equation in the form 

f(5p) = D s tep-5Mfe-5Ml + T,-5p (48) 

where 

4p = E.o.j +1 («) 

4m = E.o.j-1 (50) 

and 

T, = -D, p , - D 2 Soj-, + D 3 (; 0j - D 4 (5 1) 

The quantities £,m and Ti are constant over the iteration. 

The Newton-Raphson iteration is written in the form 

wn+,,= v n) -?| H <52) 

where n is the iteration number and where the prime denotes the derivative with respect to £,p. It 
can be shown that the derivative in the denominator takes on the same value for either > £,m or 
£p < £,m, and that it is given in finite difference form by 

f'<£p) = 2.0-D 5 |&p -S M 1-1.0 (53) 

In most cases attempted so far, the iteration converges quite rapidly, although further work is 
needed to assure that this algorithm is working properly. 
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4. Numerical Integration Computer Program Operation 

4.1. Computer Program 

The computer program ZORF4 has been written to accomplish the numerical integration 
and time-series Fourier analysis leading to a prediction of the cavity resonator impedance. 
ZORF4 is written in FORTRAN for use on the IBM-PC. This FORTRAN version could be easily 
adapted for use on other computers. 

Proper use of the program requires some insight to the methods by which the calculations 
are accomplished and how some of the input variables are interpreted and established. This 
section documents the use of the computer program. 


4.2. Program Limitations 

The dimensions of arrays in the computer program were sized based on estimates of the 
largest number of grid points in the x-direction that might be encountered for a reasonable 
problem. The acoustic particle displacement values (real, double precision) must be stored for 
only three rows of time values simultaneously. The values are written to a file on disk for each 
row as it is calculated. The maximum number of grid points in the x-direction is set at 5000 using 
the following reasoning: 

1 . Assume the maximum frequency encountered will be 40,000 Hz. 

2. The speed of sound is nominally 34,000 cm/sec (room temperature value). 

3. This gives a wavelength, X, of about 0.85 cm at 40 KHz. 

4. Assume the resonator cavity depths will be at most 10 of these shortest wavelengths (8.5 cm, 
about 3.3 in.) 

5. The minimum desired grid spacing will use at most 500 points per wavelength. 

6. Multiplying 10 wavelengths by 500 points per wavelength gives 5000 grid points, maximum 

The grid points are assumed to be equally spaced. Storage is not a problem for an IBM 
PC, so that the maximum array dimensions could be increased, if desired. 

For the FFT spectral analysis, the maximum time sequence block size is 1024 points, and 
the value chosen must be a power of 2. The time interval FFT data sampling is chosen to be an 
integral number of numerical integration time steps, say N^mp , so that every n* data point is 
stored as the time history for spectral analysis. Up to 16 frequency averages are allowed. 
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4.3. Time Constants 


4.3.1. Numerical Integration Time Step 

The numerical integration time step is related to the numerical integration stability 
condition given in Equation (16). The integration distance step size will be some small fraction of 
a wavelength required to achieve convergence. Experience in running the program indicated that 
up to several hundred points per wavelength were required to achieve acceptable convergence. 
Denote the value of points per wavelength at the shortest wavelength (highest frequency), which 
is input to the program, as Nppi. Then the distance step size is given by 


Ax = 


c 

8N PP iFph 


(54) 


where F P h is the highest frequency component in the pressure input signal and c is the speed of 
sound. The multiplier 8 in the denominator is chosen because the highest frequency in the data 
analysis (the Nyquist frequency) is set to 8 times the highest frequency in the input signal. 


The numerical integration time step, in normalized form At, is equal to Ax, by the stability 
criterion. The actual numerical integration time step will be At divided by the speed of sound, or 


At = 


1 

8N ppl F ph 


(55) 


If the minimum Ax is the shortest wavelength divided by 500, then the minimum numerical 
integration time step (in real time units) will be 


^ ~ ^ X nun _ ^-min _ ^- 8 5 

111111 ~ c “ 500c " 500*34,000 


= 5x10 8 seconds 


(56) 


where the minimum wavelength Xmm of 0.85 cm. is based on the assumption of a highest 
frequency of interest of 40,000 Hz. and a minimum speed of sound of 34,000 cm/sec. 


4.3.2. Transient Decay Time 

Before initiating data acquisition of the acoustic displacement time history for spectral 
analysis, it is necessary to be assured that the solution has achieved steady state and any transients 
caused by “switching on” the applied pressure at time zero have decayed sufficiently to have 
negligible effect. The procedure used to set the elapsed time before initiating the first data 
acquisition is described in this section. 

If we assume a maximum cavity depth of interest of 5 cm., then the longest time required 
for a plane wavefront to traverse the cavity in both directions will be 
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t _ L _ 5 cm 

^cavmax ,,. m , — 3.0E 4 sec 

c 34000 cm /sec 


( 57 ) 


If we assume that the lowest frequency of interest is 100 Hz., then the maximum wave period will 
be 


Tpcrmax = 001 Sec (58) 

These quantities can be used to estimate a required solution “settling time” during which 
transients caused by switching on the applied pressure are allowed to decay. We are assuming 
that the damping inherent in the facesheet will always be sufficient to eventually eliminate the 
transient from the solution. 

The acoustic pressure applied to the faceplate is obtained by summing a Fourier series of 
frequency components, and this pressure may not be zero at time zero. This finite jump in 
pressure at time zero will introduce a transient into the solution. To guarantee that the pressure 
starts from zero and builds slowly, a linear ramp function that extends over the first three longest 
periods (that is, the period of the lowest frequency component of the applied pressure), increasing 
from zero to one, is multiplied into the applied pressure. 

Then a settling time equal to at least three more longest periods is added, with an option 
of adding as many additional longest periods as desired. Finally, a time period equal to the time 
required for the wave to traverse the cavity (one-way) four times is added. The only means of 
determining whether, in fact, the transient has decayed acceptably is to repeat the case with an 
increased length of delay interval and compare the results with the initial run. 


4.4. Spectral Analysis 
4.4.1. Data Sampling 

When the process reaches a steady state, the data acquisition can be initiated. Time 
histories are acquired for the particle displacement and the applied pressure. Data is acquired in 
blocks in the same manner as it would be acquired by the analog-to-digital board of a digital 
spectral analyzer. FFT’s are applied to each of the data blocks and the results are then frequency- 
averaged. 

The pressure input frequencies are chosen such that they can be matched exactly by the 
data analysis output frequencies. This allows the use of a rectangular window, eliminating the 
need for more exotic windowing. 15 Some forethought is required for the initial choice of pressure 
frequencies, for they must align with the frequency bins of the FFT of the solution based on a data 
acquisition sample size that is a power of 2. 
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The numerical integration time step is much finer than the time step required for data 
sampling. The Nyquist frequency for the data analysis is set at 8 times the highest frequency in 
the pressure input signal. The sampling time interval Ax commensurate with achieving this 
Nyquist frequency is given by 


Ax 


samp 


16F, 


(59) 


ph 


The number of numerical integration time intervals per data sampling time interval is then given by 


N 


samp 


_ ^samp _ Nppl 

Ax 2 


(60) 


By choosing Nppi to be a multiple of 2, we guarantee that there will always be an integral number 
of time intervals between data sampling points. This is a convenience that makes it unnecessary 
to store the results of the numerical integration between data sampling points as long as the points 
are counted. Storing each integration point would be required if it were necessary to interpolate 
the results to the required data sampling time. 

The block size of the data acquisition, N b ik, is provided as input. This must be an integer 
power of 2 for the convenience of the FFT. The data block record length is, then, in real time 
units. 


T _ Nbik 


16F, 


ph 


and the frequency resolution of the resulting spectra is 


Af = 


16Fph 

N blk 


(61) 


(62) 


The maximum block size is Nbik = 1024, which give a finest frequency resolution of 1/64 of F ph . 
Recall that the pressure input frequencies must be input at exactly these frequencies. 

Some combinations of parameters Af, F p h , and N b ik must be chosen at the start to assure 
the data analysis frequencies and pressure frequencies match (there may be more data analysis 
frequencies than pressure frequencies, but the pressure frequencies must align with data analysis 
frequencies). If pressure data is from a measured FFT spectrum, the values may have to be 
combined in different bandwidth frequency bins to fit the data requirements of ZORF4. The 
choice of these parameters may have to be guessed iteratively, and an artificial F pb with very low 
pressure values may have to be created to obtain the desired results. 
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4.4.2. Computation of Acoustic Pressure, Velocity, and Impedance 


The acoustic pressure is needed at the surface of the perforated plate. The pressure is 
given as input to the problem, in terms of its frequency spectral components. To obtain the time 
history of the pressure at any arbitrary time, t, it is necessary to expand the Fourier components of 
the frequency, 


NFP 

p(t)=2 p F n e i<0nt (63) 

n=l 

where PF„ is the n** 1 complex spectral frequency component, ©„ is the n 111 circular frequency, and 
NFP is the number of frequencies with coefficients not equal to zero. 

The pressure on the faceplate due to a normally-incident pressure wave is assumed to be 
known. Usually, only the acoustic pressure of the incident wave would be known, but the actual 
excitation pressure is that due to the combined incident and reflected waves at the surface. The 
reflected wave component, however, depends on the impedance of the surface, which is not 
known a priori. 

If Pfwd is the amplitude of the incident (forward-traveling) pressure wave for the © 
component in the frequency domain, impinging on the facesheet from the left (outside the cavity), 
then the known forward-traveling part of the pressure wave is 

Pfwd ( x ) = P fwd el(kX_<Bt) (64) 

Using the impedance-matching boundary condition at the faceplate surface (x = 0), it can be 
shown that the resultant overall acoustic pressure at the faceplate outside surface is given by 

op 

P(0) = (65) 

1 + a 7 


where 



( 66 ) 


is the normalized acoustic admittance of the surface. Note the expected pressure-doubling when 
the impedance is infinite (hardwall surface). 

The SPL at the surface, then, must be determined from the rms magnitude of this resultant 
pressure, 
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( 67 ) 


_ ^ Pfwd nns 

rms / " 

yl + Re(a)+|a| 

The input to Z0RF4 is the complex resultant pressure p(0). 

To obtain the incident pressure wave component that will provide a desired surface SPL, 
one must work backwards iteratively through the above equations, since the impedance may 
change as a function of pressure amplitude. Given only the incident pressure wave component, 
one would have to solve iteratively to determine the impedance, since it is a function of the 
resultant pressure, which in turn depends on the impedance. 

The frequencies must be chosen to align with the frequency bins of the FFT data 
processing, and thus should be some integer multiple of Af as given in Equation (62). Note that it 
is not necessary to include a pressure spectral coefficient for each frequency between 0.0 and the 
highest frequency. The input could consist of only several widely spaced tones. The sum is over 
only those pressure coefficients with finite values, which saves some computation time. 

The input pressure coefficients PF n should be in physical pressure units of dynes/cm 2 peak 
amplitude. The SPL is then computed from 



SPL = 201 og 


p i 

V2 (0.0002) J 


( 68 ) 


where p is either the pressure in a spectral frequency bin or some combination of adjacent 
pressures. 

The spectral components of the acoustic pressure are given, but the spectral components 
of the acoustic velocity must be obtained from an FFT of the acoustic displacement time history 
solution at the faceplate. If X„ represents the n* spectral component of acoustic displacement at 
circular frequency co„, then the n^ component of acoustic velocity is given by 


V n e 


-ico.t 


= £n=-ia„X n e 




(69) 


or 


V„=-m„X n = -i2jtf 1 ,X„ 


(70) 


Given the pressure and velocity components at each frequency, the acoustic impedance 
normalized by pc is simply 


Zn PFn 

pc pcV n 


(71) 
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The real part of this is the acoustic resistance and the imaginary part is the acoustic reactance, 
both normalized by pc. 

4.5. ZORF4 Program NAMELIST Input Guide 

The computer program ZORF4 that implements the above analysis is written in 
FORTRAN code and runs on an D3M-PC compatible computer. Input to the program is through 
a NAMELIST data file, as follows: 
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ZORF4 NAMELIST INPUT GUIDE 


Time in sec, Displacements in cm., Pressure in dynes/cm A 2, Velocity in cm/sec 
DP implies a double precision variable 

&ZDAT 

HEAD ='Case Title' Filenames up to 30 characters 
DFOUT='path/filename' Case data output file 
DFTHDAT='path/filename' Time history of acquired data: 

Formatted as Block# Sample# Time Pressure Displacement(O) 
DFTRNST-path/filename' Transient time history up to start of D/A: 

Formatted as TimePnt# Time Pressure Displacment(O) 

DFSPECT- path/filename' Spectral data: 

Formatted as Freq Prms SPL Vrms Zrcn 
RLIN= Linear resistance (DC flow A-value) cgs Rayls (DP) 

RNL= Nonlinear resistance (DC flow B-value) cgs Rayls/(cm/sec) (DP) 

SIGMA= Faceplate porosity (open area ratio) (DP) 

DCAV= Cavity depth, cm. (DP) 

RHO= Air density, g/cm A 3 (DP) 

CS= Speed of sound, cm/sec (DP) 

EFFLGTH= Effective orifice length, cm (DP) 

CD= Orifice discharge coefficient (DP) 

CFNL= Nonlinear term coefficient, set to 1.D0 or 0.D0 
NFP= Number of input frequencies 

PF(1)= Complex DP array of pressures at each frequency (peak amplitude, dynes/cm A 2) 
FRQP(1)= DP array of frequencies for pressure input 
NBLKSZ= Data acquisition block size (power of 2, up to 1024) 

NTDO= Number of extra cycles for transient dieout, based on period of lowest frequency 
NDSAMP= Number of averages for frequency averaging, up to 16 

/ 
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The input data file should be given the name ZORF4IN.NML, and it should reside in the 
same directory as the ZORF4.EXE file. Be sure to back up prior NAMELIST data files that are 
to be saved with a different name. The output files are written to the directory and filename 
specified in the NAMELIST input. 


4.6. Sample Case 

A sample case was chosen to illustrate the use of the ZORF4 program. The case includes 
multiple frequency pressure input, with an applied pressure of 130 dB SPL in ten adjacent 125 
Hz. bin width frequency bands from 875 Hz. to 2000 Hz. This case will be used later to compare 
the time domain numerical integration impedance model with the frequency domain model. 

The basic input parameters for SDOF resonator are defined as follows: 

Faceplate porosity = 5.0% 

Faceplate thickness = 0.032 inches = 0.081 cm 
Faceplate hole diameter = 0. 132 inches = 0.132 cm 
Cavity depth = 1.0 inches = 2.54 cm 

Other constants are: 

Air density = 0.001206 g/cm 3 
Speed of sound = 34,394 cm/sec 
Orifice discharge coefficient = 0.76 
Absolute viscosity = 1.7894E-4 g/(cm-s) 

These are characteristic of air at about 70 degF. 

Additional inputs are the effective orifice length, the linear resistance, and the nonlinear 
resistance. These parameters, derived using equations given in the following section, are: 

Effective orifice length = 0. 176 cm 

Linear resistance = 0.5020 cgs Rayls 

Nonlinear resistance coefficient = 0.2088 cgs Rayls/(cm/sec) 

The following is an input sheet for the ZORF4 run: 

&ZDAT 

HEAD=’ SAMPLE CASE #1' 

DFOUT= 1 SCI . OUT ' 

DFTHDAT= ' ZSC1 . DAT ' 

DFTRNST= ' TSC1 . DAT ' 

DFSPECT= ' SSC1 . DAT ' 

RLIN=0. 502D0 
RNL=0.2088DO 
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SIGMA=. 05D0 
DCAV=2 . 54D0 
RHO=. 001206D0 
CS=34394 . DO 
EFFLGTH=0. 176D0 
CD=.76D0 
CFNL=1 . ODO 
NFP=10 

PF(1)=(894.43DO,O.DO) , ( 894 . 43D0, 0 . DO) , (894 . 43D0, O.DO) , (894 . 43D0, 0 . DO) 
PF(5)=(894. 43D0, O.DO) , ( 894 . 43D0, 0 . DO) , (894 . 43D0,0.D0) , ( 894 . 43D0, 0 . DO ) 
PF(9)=(894.43D0,0.D0) , ( 894 . 43D0, 0 . DO) 

FRQP ( 1 ) =875 . , 1000. , 1125. , 1250. , 1375. 

FRQP ( 6) =1500 . , 1625. ,1750. ,1875. ,2000. 

NBLKSZ=256 

NTDO=5 

NDSAMP=3 

/ 

The following is the main output run for this case, which was written into file SCI OUT : 

PROGRAM ZORF4 

TIME DOMAIN CENTRAL DIFFERENCE MODEL FOR IMPEDANCE OF 
NONLINEAR SDOF RESONATOR COVERED BY PERFORATED PLATE 

TITLE: SAMPLE CASE #1 

RLIN = .502000 cgs rayls 

RNL = .208800 cgs rayls/ (cm/sec) 

SIGMA = . 050000 

CAVITY DEPTH = 2.540000 cm 

EFFECTIVE CAVITY DEPTH = 2.536557 cm 

RHO = .001206 g/cm A 3 

CS = 34394.000000 cm/sec 

RHO*C = 41.479164 cgs Rayls 

EFFECTIVE ORIFICE LENGTH = .176000 cm 

Cd = .760000 

# OF PRESSURE FREQUENCIES = 10 

DATA ANALYSIS BLOCK SIZE = 256 

# DATA SAMPLE BLOCKS TAKEN = 3 

# INTEGRATION PNTS PER SHORTEST WAVELENGTH = 200 

# EXTRA CYCLES SETTLING TIME (NDTO) = 5 

# TIME STEPS TO D/A START = 30202 

OUTPUT DATA FILES: 

INPUT DATA AND PARAMETERS - SCI. OUT 
PRESS & VELOC TIME HIST - ZSC1.DAT 
INITIAL TRANSIENT TIME HIST - TSC1.DAT 
OUTPUT F, P, V, ZRC SPECTRA - SSC1.DAT 


J 

FREQ 

CMP LX : 

PRESSURE COEFF 

SPL 

1 

875.00 

8 . 944300E+02 

0 . 000000E+00 

130.00 

2 

1000.00 

8 . 944300E+02 

0 . OOOOOOE+OO 

130.00 

3 

1125.00 

8 . 944300E+02 

0 . 000000E+00 

130.00 

4 

1250.00 

8 . 944300E+02 

0. OOOOOOE+OO 

130.00 

5 

1375.00 

8 . 944300E+02 

0. OOOOOOE+OO 

130.00 

6 

1500.00 

8 . 944300E+02 

0. OOOOOOE+OO 

130.00 
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7 1625.00 8 . 944300E+02 0. OOOOOOE+OO 130.00 

8 1750.00 8 . 944300E+02 0. OOOOOOE+OO 130.00 

9 1875.00 8 . 944300E+02 0. OOOOOOE+OO 130.00 

10 2000.00 8 . 944300E+02 0. OOOOOOE+OO 130.00 

LOWEST PRESSURE FREQ = 875.00 Hz 

HIGHEST PRESSURE FREQ = 2000.00 Hz 

UPPER ANALYSIS FREQUENCY = 16000.00 Hz 

DATA ACQUISITION DELTA T = 3.1250D-05 sec 

TIME STEP = 3 . 12500E-07 sec 

HH or DELTA TAU = 1.07481E-02 cm 

NONLINEAR TERM COEFF = 1.0000 

CF1 = -3 . 0471D+04 

CF2 = -3 . 0470D+04 

CF3 = -6 . 0941D+04 

CF4 = 4 . 6520D+01 

CF5 = 3 . 7468D+05 

CF1 - CF4 = -3 . 0517D+04 

D1 = -2 . 2969D-11 

D2 = 9 . 9996D-01 

D3 = 1 . 9969D+00 

D4 = -3 . 0487D-03 

D5 = -1 . 2278D+01 

LONGEST CYCLE TIME = 1.1429D-03 sec 

SETTLING TIME = 9.1429D-03 sec 

CAVITY BOUNCE TIME = 7.3850D-05 sec 

TIME TO START OF D/A = 9.4381D-03 sec 


SAMPLE CASE #1 


FRQ 

PRMS 

SPL 

VHMS 

RESIS/RHOC 

REACT/ RHOC 

.000 

1. 84322E-05 

. 00 

0. OOOOOE+OO 

.00000 

. 00000 

125.000 

4 . 96309E-05 

.00 

2. 05345E-02 

- . 00006 

-.00001 

250.000 

4 . 154 85E-05 

.00 

5 . 70915E-02 

. 00000 

. 00002 

375.000 

4 . 35998E-05 

.00 

1 . 36698E-01 

. 00001 

. 00000 

500.000 

5 . 61563E-05 

. 00 

3. 11184E-01 

.00000 

. 00000 

625.000 

7 . 44 8 67E-05 

. 00 

6. 98051E-01 

. 00000 

.00000 

750.000 

1 . 322 8 0E-04 

.00 

1 . 51042E+00 

.00000 

.00000 

875.000 

6 . 32457E+02 

130.00 

8 . 00341E+00 

.77985 

-1.73821 

1000.000 

6 . 32457E+02 

130.00 

9 . 62628E+00 

.85589 

-1.33280 

1125.000 

6 . 32458E+02 

130.00 

1. 13897E+01 

.88111 

-1.00787 

1250.000 

6 . 3245 8E+02 

130.00 

1. 32377E+01 

.88617 

-.73581 

1375.000 

6 . 32457E+02 

130.00 

1 . 53689E+01 

.85339 

-.50597 

1500.000 

6 . 32457E+02 

130.00 

1. 82088E+01 

.78394 

-.29433 

1625.000 

6 . 32457E+02 

130.00 

2 . 36269E+01 

.64530 

-.00767 

1750.000 

6 . 32458E+02 

130.00 

2. 07269E+01 

.67415 

.29444 

1875.000 

6 . 32458E+02 

130.00 

1. 89357E+01 

.64474 

.48240 

2000.000 

6 . 32458E+02 

130.00 

1 . 75159E+01 

.54415 

.67946 

2125.000 

2. 93880E-04 

3.34 

7 . 70186E+00 

.00000 

.00000 


The results of the computation are found in the last table. A plot of the calculated impedance is 
shown in Figure 6. 
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The convergence of the numerical integration solution depends critically upon the number 
of grid points per wavelength used for the computation. Table 1 below shows the solution 
convergence comparing calculations with 100 and 200 points per wavelength. Notice that the 
solution has converged out to the third decimal place. 



100 pts/wavelength 

200 pts/wavelength | 

Frequency 

Resistance 

Reactance 

Resistance 

Reactance 

875 

0.7811 

-1.72875 

0.77985 

-1.73821 

1000 

0.85565 

-1.32442 

0.85589 

-1.3328 

1125 

0.88038 

-0.9995 

0.88111 

-1.00787 

1250 | 

0.88628 

-0.72759 

0.88617 

-0.73581 

1375 

0.85283 

-0.49994 

0.85339 

-0.50597 

1500 

0.77984 

-0.28886 

0.78394 

-0.29433 

1625 

0.64522 

0.00498 

0.6453 

-0.00767 

1750 

0.6785 

0.30041 

0.67415 

0.29444 

1875 

0.64803 

0.48742 

0.64474 

0.4824 

2000 

0.54877 

0.68181 

0.54415 

0.67946 


Table 1. Comparison of ZORF4 impedance computation convergence for 100 and 200 
points per wavelength. 
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Z0RF4 also creates four auxiliary data files. The most useful file is the ASCII table of 
pressure, velocity, and impedance versus frequency, for it can be imported easily into a plotting 
program. Users should be aware that the transient time history data files can become several 
megabytes large. Unless the transient behavior is of interest, these files should be deleted. 


4.7. Check for Transient Decay 

Assurance is needed that transients in the solution caused by “switching on” the pressure 
excitation at time zero have decayed to the point of having negligible effect on the impedance 
determination. It is assumed that the natural damping in the facesheet will act to attenuate these 
transients in the numerical integration. 

The time constant used to attempt to achieve transient elimination is the period of the 
lowest frequency component in the pressure excitation (assuming no zero-frequency component). 
The pressure excitation “ramps up” linearly over the first three (or more) longest periods plus four 
cavity traverse times as the default settling time condition. 

The user has an option of adding an arbitrary number of additional long period cycles to 
the settling time. The sample case defined in Section 4.6 was run with two different transient 
settling times, first with no extra periods in the settling time (default condition), and then with an 
additional five longest periods added, and the results were compared. 

Table 2 shows the comparison of the impedance calculation at each frequency for the two 
different transient decay times. The excellent agreement can be noted, indicating acceptable 
convergence for this case. The actual time history of the acoustic particle displacement transient 
for the default settling time case, up to the time of initiation of data acquisition, is shown in Figure 

7. 



Normalized Impedance | 

Frequency, Hz. 

Decay time = 0.00372 sec 

Decay time = 0.00944 | 

875 

1.0593 

-1 .7051 i 

1.0629 

-1.7007i 

1000 

1.1487 

-1.2332i 

1.1491 

-1.2396i 

1125 

1.1681 

-0.8660i 

1.1652 

-0.861 3i | 

1250 

1.1296 

-0.5240i 

1.1365 

-0.5275i 

1375 

1.0632 

-0.1864i 

1.0511 

-0.1 842i 

1500 

1.0705 

0.2379i 

1.0862 

0.2333i 

1625 

1.1525 

0.4354i 

1.1416 

0.4420i 

1750 

1.0864 

0.651 1i 

1.0930 

0.6485i 

1875 

0.9971 

0.8686i 

0.9916 

0.8667i 

2000 

0.8238 

1 . 1 0521 

0.8258 

1.1094S 


Table 2 Comparison of impedance computation for two transient decay times. 
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Figure 7. Transient of acoustic particle displacement for sample case, default settling 
time. 

Since the Fourier transforms and averaging were taken for two different samples, some of 
the difference may be due to differences in the time samples, as opposed to being caused by the 
transient. In any case, neither cause has led to large differences. Further parametric studies are 
required before concluding unequivocally that there is no transient decay problem or to develop 
generalized acceptable decay time rules. 
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5. Comparison of Time-Domain Prediction to Frequency-Domain Models 

5.1. Relation of Parameters to Current Impedance Models 

The quantities m, Rl, and Rnl are related to the physical parameters in the perforated 
facesheet. Conventionally-used models are discussed in detail in Reference 16 . In this section, we 
will first consider the mass reactance and then the resistance. 


The mass reactance of a perforated sheet is given by 

X m _ 27cf(t + sd) 
pc caC d 


(72) 


where t is the facesheet thickness, d is the hole diameter, and e is an end correction coefficient. 
The effective orifice length is given by 


f = t + ed (73) 

Various expressions for the end correction coefficient e are discussed in Reference 16. 

From Equation (28), we can identity m as a factor in the mass reactance, such that 

v 2nf , , i 

X m = m = -km = k (74) 

c aC d 


It is important to note that the expression for mass reactance in Equation (74) differs from that in 
Reference 16 by the presence of the orifice discharge coefficient in the denominator. Most 
derivations of mass reactance for a perforate omit this factor. Since the effective orifice length is 
usually empirically-determined, any effect of the 1/C d factor is incorporated into the effective 
orifice length. 

The mass reactance is known to vary with the acoustic velocity amplitude (that is, with 
SPL) as well as being proportional to frequency. The variation with porosity may be more 
complex than a simple inverse variation. In order to be able to model the effects on mass 
reactance of the velocity-dependence of the slug of air in the orifice (the effective length) and the 
possibly independent variation of the orifice discharge coefficient separately, the form that has C d 
in the denominator is chosen here. 


A commonly used expression for the linear part of the DC flow resistance is 


r l = 


32pot 

oC d d 2 


(75) 
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where v is the kinematic viscosity. A model for the coefficient of the nonlinear part of the DC 
flow resistance is 


R-nl ~ 


2(0 C D )' 


(76) 


Some interpretations are necessary to adapt these models for use in Z0RF4. The 
empirical rule for the equivalence of DC flow resistance measurements and acoustic faceplate 
resistance is to identify the DC flow incident velocity with the rms acoustic incident velocity. 

From Equation (1), the DC flow resistance across the facesheet is given by 


r =^Edc =A + BU 

DC JJ 


(77) 


where U is the incident DC flow velocity incident upon the facesheet. The coefficients A and B 
are the measured linear and nonlinear DC flow resistance coefficients, respectively. The 
corresponding acoustic resistance is given by 


Rag - R LIN + R NL u nns ( 78 ) 

where Unns is the rms acoustic particle velocity. Assuming the equivalence of the resistance for the 
DC flow and the acoustic cases, and the equivalence of U and u™ , we identify Run = A and 
Rnl = B. 

The frequency-domain impedance models use the assumption that Unns can be identified 
with the DC flow velocity to determine the faceplate resistance at any given frequency. Due to 
what is called a “bias effect”, it has been found that the validity of the frequency-domain models 
can be improved by makin g the additional assumption that u^ is the overall rms acoustic velocity, 
integrated over the entire frequency spectrum. This leads to a predicted resistance that is 
essentially constant with frequency, since the overall Unns is the same at each frequency. 

The assumption for the time-domain model is that the instantaneous flow resistance across 
the facesheet is given by the DC flow resistance relationship where U is identified with the 
instantaneous value of velocity as computed by the numerical integration. In this manner, the 
faceplate velocity contains all frequency components. This may reproduce the frequency-domain 
bias effect, but should maintain any effects of variation in pressure frequency spectrum shape. 
The issue is whether the DC flow measurement coefficients A and B are appropriate for the time- 
varying acoustic flow, or whether more sophisticated models are required. This issue must 
ultimately be resolved by comparison with experimental data. 

The objectives for developing the time-domain model are to confirm or improve the 
perforated faceplate resistance and mass reactance models based on insight gained from the 
results of the numerical computations and their comparison to measured data. The numerical 
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model allows fine-tuning of the effects of variation of any of the parameters that make up the 
impedance model, such as the orifice discharge coefficient or effective orifice length. 


5.2. Comparison of Time-Domain and Frequency-Domain Predictions 

5.2.1. Multiple-Frequency Nonlinear Sample Case 

The impedance of the treatment panel defined in the sample case in the previous section 
was predicted using a current frequency-domain model. Figure 8 shows a comparison of the 
impedance predicted by the time-domain model and the frequency-do main model. 



Note that, for this particular case, the time-domain model predicts both higher resistance 
and reactance. The added reactance is caused primarily by a higher prediction of mass reactance 
by the time-domain model. This reflects the use of the Ca factor in the deno min ator of the time- 
domain model and its absence in the frequency-domain model. Note that the resistance predicted 
by the time-domain model is relatively flat, but not perfectly flat like the prediction for the 
frequency-domain model. 

As we shall see later, one cannot at this point make the generalization that the time- 
domain model always over-predicts the resistance or mass reactance. Each case must be 
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examined individually. The ultimate resolution of model accuracy must come from comparison 
with experimentally-measured values, and this will be the subject of the following section. 


5.2.2. Linearized Multiple-Frequency Sample Case 

In the previous section, we found a disagreement between the time-domain and frequency- 
domain impedance models both with regard to nonlinear resistance and mass reactance. 
Assurance is needed that the two models give equivalent results for a case where the 
discrepancies caused by the nonlinear resistance iteration and the assumptions about the presence 
of the discharge coefficient in the mass reactance denominator are eliminated. 

In this section, a sample case is considered in which the resistance is purely linear and the 
discharge coefficient is included in the mass reactance term for both models. The faceplate 
parameters remain the same as the previous sample case. 

To linearize the resistance in the frequency-domain model, the DC flow resistance B- 
coefficient is set arbitrarily to zero and a grazing flow effect for Mach 0. 1 is added to give a total 
resistance of 1 .017 pc. In the time domain model, then, the Run input is set equal to 1.017 pc and 
the Rnl input is set equal to zero. The results of the predictions are shown in Table 3 . 



Frequency Domain 

Time Domain 

Frequency 

Resistance 

Reactance 

Resistance 

Reactance 

875 

1.0169 

-1.5859 

1.0199 

-1.5782 

1000 

1.0169 

-1.1522 

1.0202 

-1.1455 

1125 

1.0169 

-0.7866 

1.0199 

-0.7805 

1250 

1.0169 

-0.4688 

1.0202 

-0.4631 

1375 

1.0169 

-0.1854 

1.0199 

-0.1802 

1500 

1.0169 

0.0721 

1.0201 

0.0773 

1625 

1.0169 

0.3099 

1.0200 

0.3145 

1750 

1.0169 

0.5324 

1.0199 

0.5372 

1875 

1.0169 

0.7426 

1.0201 

0.7469 

2000 

1.0169 

0.9431 

1.0198 

0.9477 


Table 3. Comparison of linear impedance predicted by frequency-domain and time- 
domain methods. 

Note the excellent agreement, assuring that there are no fundamental differences between the 
frequency-domain computation and the time-domain computation, at least for this linear sample 
case. 
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6. Comparison of Predicted and Measured Impedance 

6.1. Single-Frequency Pressure Excitation 

Normal incidence impedance tube measurement data were obtained for two perforated 
plate facesheets using pure tone excitation of a known SPL. The two perforates differ primarily 
in their porosity, and are defined in Table 4. 


Parameter 

Faceplate 1 

Faceplate 2 

Porosity 

0.045 

0.157 

Thickness, inches 

0.032 

0.032 

Hole diameter, inches 

0.060 

0.050 

Cavity depth, inches 

0.75 

0.75 

Frequency, Hz. 

1500 

3000 

Linear Resistance, cgs Rayls 

0.594 

0.245 

Nonlinear Resistance, 
cgs Rayls/cm 

0.526 

0.0432 

Effective Orifice Length, cm 

0.192 

0.159 


Table 4. Facesheet definition for single frequency predicted and measured impedance 
comparison. 

Two SPL values are used in each case, 130 dB and 140 dB. Tables 5 and 6 are 
comparisons of impedances predicted using the frequency-domain model both with and without 
the Cd factor in the denominator of the mass reactance, the time-domain model, and the measured 
value. 


Plate 1 - 4.5% Porosity 

Normalized Impedance I 

Model 

130 dB 

140 dB 

Measurement 

0.587 - 0.222i 

1.874 -0.489i 

Time-domain model 
with C d in mass reactance 

0.464- 0.1 77i 

0.839 - 0.165i 

Frequency-domain 
with C d in mass reactance 

0.416 -0.201i 

0.761 - 0.201 i 

Frequency-domain 
no C d in mass reactance 

0.303 - 0.570i 

0.678 - 0.570i 


Table 5. Comparison of predicted and measured impedance for Faceplate #1. 


Plate 2 - 15.7% Porosity 

Normalized Impedance ! 

Model 

130 dB 

140 dB 

Measurement 

0.188 + 0.1248i 

0.330 + 0.069i 

Time-domain model 
with C d in mass reactance 

0.104 + 0.1 55i 

0.221 + 0.1 55i 

Frequency-domain 
with C d in mass reactance 

0.092 + 0.1 51 i 

0.199 + 0.151i 

Frequency-domain 
no C d in mass reactance 

0.123 -0.024i 

0.222 - 0.024i 


Table 6. Comparison of predicted and measured impedance for Faceplate #2. 
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The comparisons show mixed results. Generally, all models tend to under-predict the 
measured resistance. For both cases, the reactance values tend to bracket the measured value. 
No clear advantage can be seen for either the time-domain or the frequency-domain models, but it 
is also obvious that the time-domain model is not generating large errors. One can conclude that 
both models require improvement in the modeling of the resistance and mass reactance of the 
perforate if increased accuracy is the objective. 


6.2. Multiple-Frequency Pressure Excitation 

An acoustic impedance measurement was performed by Rohr, Inc using a perforate plate 
SDOF treatment with the following parameters: 

Porosity = 4.43% 

Hole diameter = 0.0453 inches = 0.1151 cm. 

Faceplate thickness = 0.03 1 inches = .07874 cm. 

Cavity depth = 0.75 inches = 1.905 cm. 

The measurement temperature was 67.9 degF. The treatment was excited by a fairly flat 
broadband pressure spectrum with an overall SPL of 149 dB. The time-domain prediction was 
made by combining the original measured 8 Hz. binwidth measured spectral pressure values into 
far fewer 576 Hz. bins. This was required by numerical computation problems encountered in the 
ZORF4 runs. 

Figure 9 is a comparison of the measured impedance with the predictions from the 
frequency-domain and the time-domain models. The following observations can be made: 

• The frequency-domain model predicts the measured resistance quite closely. The measured 
impedance is quite flat, although some variation not predicted by the frequency-domain model 
can be observed. 

• The frequency-domain model overpredicts the mass reactance by a significant amount. This is 
in part due to the inclusion of the 1/Cd term in the mass reactance for this prediction. This 
factor was included to be able to compare the results directly to the time-domain prediction. 

• The time-domain model overpredicts the measured resistance by a large amount above about 
1700 Hz. 

• The time-domain model shows a large overprediction of the mass reactance. 
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Figure 9. Comparison of measured and predicted frequencies for SDOF treatment 
with broadband pressure excitation. 

The extreme overprediction of resistance and mass reactance by the time-domain model 
appears to be characteristic of the computation for high SPL values and multiple frequency input. 
The problem is exacerbated by large numbers of pressure frequencies, under which condition the 
computation will go unstable. The problem is caused by a combination of the nonlinear resistance 
algorithm, the high SPL levels, and the large number of pressure frequencies, because it is not 
experienced for the same case but with a linearized resistance. Further work is needed to examine 
the behavior of the nonlinear resistance algorithm. 


38 



7 . 


Recommendations for Further Research 


A more efficient finite difference integration algorithm is needed both to improve accuracy 
and reduce computation time. The objective of the development of this program, to obtain an 
operational computation under the simplest of analytical conditions, has been met, but this 
upgrade is needed to improve the utility of the program. Numerous improved algorithms are 
available, and one that should receive strong consideration is the dispersion-relation-preserving 
finite difference scheme of Tam 17 ’ 18 . 

The validity of the algorithm for the nonlinear impedance effects needs to be checked 
further, particularly for multiple-frequency cases. Possibly an improved finite difference algorithm 
may improve the nonlinear iteration, as well. This effort is necessary before the program can be 
used to examine effects of pressure frequency spectra shape. 

Rather than expand the pressure Fourier series at each numerical integration time point, it 
may be possible to use Inverse Fast Fourier Transform techniques to generate the pressure time 
history for a complete data acquisition block. If this is possible with sufficient accuracy, this 
would speed the computation significantly. 

The development and incorporation of more advanced resistance and mass reactance 
models for incorporation into the model is needed. In particular, the Reynold’s number 
dependence of the orifice discharge coefficient and the velocity dependence (nonlinearity) of the 
mass reactance could be incorporated. 
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8 . 


Nomenclature 


A 

B 

c 

C d 

Ci-C 6 

d 

D 

D,-D 5 

F* 

h 

i 

j 

k 

l 

L 

m 

M 

N 

Nbik 

NFP 

Nfr, 

Nppi 

N 

■f^samp 

P 

PF„ 

Pfwd 

Pfwd 

R 

Rl 

Rnl 

SPL 

t 

T, 

Tree 

u 

Hons 

V 

Vi 

V n 

X 

X„ 


linear DC flow resistance coefficient 
nonlinear DC flow resistance coefficient 
speed of sound 
orifice discharge coefficient 

numerical integration difference equation coefficients 
perforate hole diameter 

effective acoustic damping per unit area of faceplate 

numerical integration difference equation coefficients 

highest frequency present in pressure excitation spectrum 

numerical integration time step 

square root of -1 or integer index 

integer index 

wave number 

effective orifice length 

cavity depth 

mass reactance coefficient 

effective acoustic mass per unit area of faceplate 

number of spatial grid points in cavity 

block size for data acquisition 

number of non-zero pressure frequencies 

number of frequency bins in frequency spectrum 

number of points per shortest wavelength in x-grid 

number of numerical integration time intervals between each data sample 

acoustic pressure 

spectral component of applied pressure at n* frequency 

forward-traveling incident pressure wave amplitude 

forward-traveling part of incident pressure wave 

resistance coefficient 

linear part of facesheet damping 

nonlinear part of facesheet damping 

Sound Pressure Level, dB, re: 0.0002 dynes/cm 2 

Facesheet thickness, also time 

interim variable for nonlinear numerical integration iteration 
data block record time length 

incident flow velocity for DC flow resistance measurement 

rms value of acoustic velocity 

acoustic particle velocity 

acoustic velocity in the i* frequency band 

n 1 * 1 spectral component of acoustic velocity 

spatial variable along cavity 

n* spectral component of acoustic displacement 

n* spectral component of acoustic impedance 
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a normalized acoustic admittance 

s mass reactance end correction coefficient 

Af Fourier analysis frequency resolution, bin size 

At grid time step size 

Ax grid spatial step size 

X wavelength 

p absolute viscosity 

v kinematic viscosity, p/p 

p medium density 

ct perforate open area ratio, or porosity 

t numerical integration time step parameter, = ct 

© circular frequency 

acoustic particle displacement 

interim variable for nonlinear numerical integration iteration 
interim variable for nonlinear numerical integration iteration 


41 



References 


Motsinger, R. E., and Kraft, R. E., “Design and Performance of Duct Acoustic 
Treatment”, Chapter 14 in Hubbard, Harvey H., ed., Aeroacoustics of Flight Vehicles: 
Theory and Practice, Volume 2: Noise Control , NASA Reference Publication 1258 Vol 
2, 1991. 

Rice, Edward J., “A Model for the Acoustic Impedance of a Perforated Plate Liner with 
Multiple Frequency Excitation”, NASA TM X-67950, October, 1971. 

Cummings, A., “Transient and Multiple Frequency Sound Transmission Through 
Perforated Plates at High Amplitude”, AIAA 84-23 1 1 , October, 1984. 

Zorumski, William E. and Parrott, Tony L., “Nonlinear Acoustic Theory for Thin Porous 
Sheets”, NASA SP-189, October, 1968. 

Zorumski, William E. and Tester, Brian J., “Prediction of the Acoustic Impedance of Duct 
Liners”, NASA TM X-73951, September, 1976. 

Rice, Edward J., A Model for the Pressure Excitation Spectrum and Acoustic Impedance 
of Sound Absorbers in the Presence of Grazing Flow”, AIAA 73-995, October, 1973. 
Hersh, A. S. and Rogers, T., “Fluid Mechanical Model of the Acoustic Impedance of 
Small Orifices”, AIAA 75-496, March, 1975. 

Hersh, Alan. S. and Walker, Bruce, “Fluid Mechanical Model of the Helmholtz 
Resonator”, NASACR-2904, September, 1977. 

Hersh, A S., “Nonlinear Behavior of Helmholtz Resonators”, AIAA 90-4020 October 
1990. 

Cummings, A. and Eversman, W. “High Amplitude Acoustic Transmission Through Duct 
Terminations: Theory”, Journal of Sound and Vibration, 91, 1983, pp.503-518. 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P., Numerical Recipes 
in FORTRAN. The Art of Scientific Computing, Cambridge University Press, Second 
Edition, 1992, pp. 825-831. 

Berezin, I. S. and Zhidkov, N. P., Computing Methods: Volume II, Pergamon Press, 1965, 
pp. 470-486. 

Smith, G. D., Numerical Solution of Partial Differential Equations, Oxford University 
Press, 1965, pp. 111-114. 

Fox, L., Numerical Solution of Ordinary and Partial Differential Equations, Pergamon 
Press, 1962, pp. 218-222. 

Lynn, P A. and Fuerst, W., Digital Signal Processing with Computer Applications, John 
Wiley and Sons, 1989, p. 251. 

Motsinger, R. E., and Kraft, R. E., “Design and Performance of Duct Acoustic 
Treatment”, Chapter 14 in Hubbard, Harvey H., ed., Aeroacoustics of Flight Vehicles: 
Theory and Practice, Volume 2: Noise Control, NASA Reference Publication 1258 Vol 
2, 1991, pp. 180-181. 

Tam, Christopher K. W., and Webb, J. C., “Dispersion-Relation-Preserving F ini te 
Difference Schemes for Computational Acoustics”, Journal of Computational Physics 
107, 1993, pp. 262-281. 

Tam, C. K. W., “Computational Aeroacoustics: Issues and Methods”, AIAA 95-0677, 
January, 1995. 



REPORT DOCUMENTATION PAGE 

Form Approved 

QMB No. Q 7Q4-Q188 

Public reporting burden for this collection of information is estimated to average t hour per response, including the time for reviewing instructions, searching existing data 
sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other 
aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and 
Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), 
Washington, DC 20503. 

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 

April 1999 Contractor Report 

4. TITLE AND SUBTITLE 

Acoustic Treatment Design Scaling Methods 

Volume 4: Numerical Simulation of the Nonlinear Acoustic Impedance of a Perforated 
Plate Single-Degree- of-Freedom Resonator Using a Time-Domain Finite Difference 
Method 

5. FUNDING NUMBERS 

C-NAS3-26617 
TA 25 

WU-538-03- 12-02 

6. AUTHOR(S) 

R. E. Kraft 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

General Electric Aircraft Engines (GEAE) 

P.O.Box 156301 
Cincinnati, OH 45215-6301 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, V A 23681-2199 

10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

N AS A/CR- 1 999-209 1 20/VOL4 

1 1 . SUPPLEMENTARY NOTES 

Lewis Project Manager: Christopher E. Hughes 

Langley Technical Monitors: Tony L. Parrott, Lorenzo R. Clark 

Prepared for Langley Research Center under Contract NAS3-26617, Task 25. 

12a. DISTRIBUTION/A VAILABILITY STATEMENT 

Unclassified-Unlimited 

Subject Category 71 Distribution: Standard 

Availability: NASA CASI (301) 621-0390 

12b. DISTRIBUTION CODE 


13. ABSTRACT (Maximum 200 words) 

Single-degree-of- freedom resonators consisting of honeycomb cells covered by perforated facesheets are widely 
used as acoustic noise suppression liners in aircraft engine ducts. The acoustic resistance and mass reactance of 
such liners are known to vary with the intensity of the sound incident upon the panel. Since the pressure drop 
across a perforated liner facesheet increases quadratically with the flow velocity through the facesheet, this is 
known as the nonlinear resistance effect. In the past, two different empirical frequency domain models have 
been used to predict the Sound Pressure Level effect of the incident wave on the perforated liner impedance, one 
that uses the incident particle velocity in isolated narrowbands, and one that models the particle velocity as the 
overall velocity. In the absence of grazing flow, neither frequency domain model is entirely accurate in 
predicting the nonlinear effect that is measured for typical perforated sheets. The time domain model is 
developed in an attempt to understand and improve the model for the effect of spectral shape and amplitude of 
multi -frequency incident sound pressure on the liner impedance. A computer code for the time-domain finite 
difference model is developed and predictions using the models are compared to current frequency -domain 
models. 


14. SUBJECT TERMS 

Aircraft noise; acoustic treatment; fan noise suppression; scale models; 
acoustic impedance 

15. NUMBER OF PAGES 

50 

16. PRICE CODE 

A03 

17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 

18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 

19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 

20. LIMITATION “ 

OF ABSTRACT 

UL 


N5N 7540-01 -280-5500 Standard Form ^b (Rev. 2-89) 

Prescribed by ANSI Sta. Z-39-18 
298-102 

















